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In a recent Letter we introduced Hellmann-Feynman operator sampling in diffusion Monte Carlo 
calculations. Here we derive, by evaluating the second derivative of the total energy, an efficient 
method for the calculation of the static density-response function of a many-electron system. Our 
analysis of the effect of the nodes suggests that correlation is described correctly and we find that 
the effect of the nodes can be dealt with. 
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I. INTRODUCTION 

Diffusion Monte Carlo (DMC) represents a powerful 
method for the accurate computation of properties of 
molecules and solids However, so far few attempts 
0, H| have been made to use DMC to calculate the static 
density-response function [4j, which is a central quan- 
tity in the analysis of many-electron systems and time- 
dependent density-functional theory [5[. One reason is 
the technical difficulty inherent in the most straightfor- 
ward method to do so: For a given perturbing potential 
one calculates the total energy at different strengths and 
numerically determines the second derivative. This then 
gives a DMC estimate of the diagonal term of the static 
response function x- There are, however, several obvi- 
ous difficulties with this. One needs one loop for various 
perturbation strengths, another loop for each k that one 
wishes to sample, and if one wants the off-diagonal terms 
a third loop for the kl . Inside each of these loops sits 
an entire wavefunction reoptimization cycle and a com- 
plete DMC run. The perturbations must be small enough 
not to change the wavefunction qualitatively and large 
enough to allow for sensible numerical derivatives. 

In a recent Letter JfJ], we showed how "applying" the 
Hellmann-Feynman 7] (HF) derivative to the DMC al- 
gorithm leads to a new algorithm, Hellmann-Feynman 
sampling (HFS), that correctly samples the first deriva- 
tive of the energy, i.e., an expectation value of an oper- 
ator. HFS works because DMC yields the correct total 
energy for nodes defined by the trial wavefunction. For 
technical reasons the operators sampled must be diago- 
nal in real space. Extending the analysis to the second 
derivative yields a DMC algorithm for the fixed-node (fn) 
static density-response function. Note that even for a 
trial wavefunction with correct nodes the fixed-node den- 
sity response is not the exact value as the real response 
includes effects from the change of the nodes. However, 
comparison with Ref. Q where the nodal variation of an 
underlying Kohn-Sham (KS) system is implicitly used, 
shows that these effect can be accounted for by general- 



izing the RPA analysis [8] to fn systems. The resulting 
method can be performed within a single DMC run, and 
in the case of inhomogeneous systems can produce off- 
diagonal elements of x as easily as the diagonal terms. 

The present paper is organized as follows. After a 
brief recapitulation of HF sampling, we derive formulae 
for the DMC sampling of x along the same line. We 
then briefly discuss technical aspects (convergence with 
respect to population size, time step, etc.). Finally, we 
look at the density response of the interacting and non- 
interacting uniform electron gas, analyze the effects of 
the nodes, and compare our results with the literature. 
Our method should also enable DMC calculations of the 
static-response function of real solids, never done before. 
We use atomic units throughout. 



II. HELLMANN-FEYNMAN SAMPLING AND 
THE DENSITY RESPONSE FUNCTION 

A. Application to the second derivative of the 
energy 

Fixed-node DMC yields by construction 
the normalized expectation value (O)dmc = 



(^ T \0\% n )/(^ T \% n ), where is the ground- 

state wavefunction constrained by the nodes of the 
Fermionic many-body trial wavefunction vf^; HFS 

correctly calculates (*o™|0|*o"}/(*o"l*o") while 
maintaining the basic DMC algorithm that samples 
^t^o"- This is because the total energy is evaluated 
correctly within standard DMC, and crucially operator 
expectation values can be cast as HF derivatives of 
the total energy. Keeping in mind that ultimately the 
DMC algorithm is nothing but a large sum that yields 
the total energy, we see that the HF derivative can be 
applied to the algorithm itself! One advantage over 
numerical derivatives is that the resulting formula can 
handle several operators simultaneously in a single DMC 
run, and maintaining orbital occupancy for perturbed 
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Hamiltonians ceases to be a problem. The DMC 
algorithm only involves numbers, so non-commutability 
of operators is no problem. Writing down the DMC 
algorithm as a mathematical formula and applying the 
HF derivative to it yields an object that when sampled 
using standard DMC produces the exact operator 
expectation value. We find that given a Hamiltonian 
H(a) — H + aO, evaluating the growth estimator of the 
energy E GR at time step i to first order in a yields a 
growth estimator that samples the operator O. Similarly 
the direct estimator E of the energy yields another 
estimator. These are Eqs. (8) and (9) of Ref. 0: 



Of R = 
Of = 
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(1) 
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da 



Of - t (EfXi - El ■ Xi) (2) 



where the bar refers to the DMC sampling at time step 



»: X = £ 



N„ 



is the total weight of walker 



h X.j = iELiOfcj. and °k,j is ^ evaluated for 
walker j at time step k. Now we assume two perturba- 
tions of the form a&A and (30b, and following Ref. @ 
we obtain growth and direct estimators of the response 
function xab — ^§j$- = from the second derivative 
of the growth and direct estimators of the energy: 



-t [xfx? -xf-x?] 

-t[XfOf-Xf-Of- 



OfX 



(3) 

of-xf] 



-?[{Ef ~ Ef){Xf - Xf){Xf - X?)] . (4) 



Any number of operators O a and O b can be sampled in 
parallel within a single DMC run. From now on we use 
the Fourier components of the density sinfca; and cosfca;. 
Note also that,as in the case of the first derivative dis- 
cussed in Ref. [f| , the growth estimator at i is already an 
averaged quantity. This property makes it an attractive 
choice for a DMC calculation. Equation Q is not only a 
more complicated formula than Eq. @ , it also has to be 
summed at each time step i. If we wish to sample many 
components of x at the same time, a large array with a 
size quadratic in the number of components of x m Eq. 
(UJ has to be generated and dealt with at each time step. 
By contrast, Eq. §3§ only has to be built at whatever 
time step one wishes to calculate x- This means that at 
each time step one only has to maintain the xf^ B , which 
is less memory intensive and much faster to compute. 
Hence we shall only use the growth estimator Eq. 



B. Computational implementation 

The growth estimator has the advantage that for each 
time step and walker we only need to deal with simple 
sampling of Nk variables for the components of the den- 
sity. The entire density-response function, including off- 



diagonal elements, can then be calculated as a correla- 
tion function Q of these variables at the end of the run 
saving computer time and memory. As in HFS, we find 
that noise rises as the sampling progresses, thus limiting 
the statistical error of the final result even if the sam- 
pling is continued indefinitely. So to reduce statistical 
noise, we increase the number of walkers instead. This 
has the additional advantage of reducing any population 
bias. We converged this by looking at population sizes of 
200, 1000, 5000, and 50000 walkers. We used the latter 
for the main results shown here (at N = 114 electrons). 
Looking at different time steps, we found that too large a 
time step shows up as a levelling off of x at a finite value 
at larger k instead of showing the correct 1/k 2 behavior. 
Interestingly, even at 8t = 0.1 the calculated x remained 
unbiased up to and well beyond k = 5kp. A large time 
step is desirable as equilibration will be faster. Here we 
use 5t = 0.01. To monitor equilibration we artificially 
extract a value Xi that when summed over all time steps 
gives the growth estimator x(iV) at the time t = NAt of 
sampling: From 



w-i 



N 



and 



1 —j2m=x(N-i) 



1 N 



it follows that 



X (N) = N XN - (N - 1) X (N - 1) 



(5) 



(6) 



(7) 



Following, e.g. three typical k's as the sampling pro- 
gresses we found that x converges exponentially. A quick 
run using few walkers in a smaller system can then be 
used to roughly estimate the convergence time which one 
then uses in an actual run. Convergence can be improved 
by setting up the sampling such that one ignores the im- 
plied x during equilibration. In order to do that it is not 
necessary to actually reverse engineer these x f° r each 
element of the density response function. Once one has 
decided on an equilibration time, N eq the desired result 
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N -N, 



eq . 



i 



i=N c „ + l 



N — N, 



[N X (N) - N eqX (N eq )} 
(8) 

It is thus sufficient to perform the costly sampling of the 
full growth estimator only twice during the run, once af- 
ter equilibration, and once at the end of the run. Doing 
so, efficiently removes much of a 1/A-like term without 
incurring much extra computation. We found that equi- 
libration was essentially independent of St and seemed 
to depend only weakly on the population size. We note 
that while in this paper we only calculate the diagonal 
terms of x> sampling the full density-response function, 
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FIG. 1: The dashed line shows the density response function 
Xo of an infinitely large unpolarized non-interacting homoge- 
neous electron gas (Lindhard function) at r 3 = 2. The black 
lines close to the Lindhard function show the exact finite-size 
Lindhard function Xo S with 114 and 4218 electrons, the latter 
following the Lindhard function closer. The dotted line shows 
the RPA response function \rpa, and the three remaining 
dot-dashed lines show (from top to bottom) the fixed-node 
Lindhard function Xo" a ^ H4 electrons, the corresponding 
fixed-node density response of an interacting system x^ n > an d 
the fixed- node RPA Xrpa- The "wiggles" are not noise as 
they correspond to the shell structure seen in the exact non- 
interacting finite-size Xq S 



including all the off-diagonal terms, is not more difficult: 
All that is needed is the evaluation of the correlation 
function of XkS at different k and k 1 . 



III. RESULTS 
A. The system 

In order to demonstrate our method we calculated the 
diagonal terms of the static density response function of 
an unpolarized free electron gas for electron-density pa- 
rameter r s = 2, 5, and 10 and corresponding density no- 
We set up the DMC calculation using 114 electrons in a 
simple-cubic super cell. We also looked at fee unit cells 
and smaller systems with 66 electrons, however, we found 
no qualitative difference. In all, we calculated the den- 
sity response at all 119 independent k- vectors between 
k = and k = 5kp. Our DMC calculations employed 
trial wave functions "Ft of the Slater- Jastrow type with 
a standard correlation term. Prior to the DMC run >Ft 
was optimized in a variance minimization run. We used 
the CASINO [3 code for all our computations. 



B. The density-response function 

In general, the response function is given by 

,, B = f2 S (°'°f<'g*'°> , (9) 

i=l 1 

where the sum runs over all excited states of the many- 
electron system, fn DMC yields the ground-state energy 
for nodes given by the trial wave function. Therefore the 
second derivative yields the fixed-node response x^™ °f 
a system for which the nodes are the same for all per- 
turbing potentials. Since in the case of a fixed-node sys- 
tem the sum entering Eq. (J9j) runs over a set of fixed- 
node excited states that differ from the actual excited 
states of the many-electron system, the fixed-node and 
non fixed-node non-interacting density-response function 
differ considerably (Fig.[lJ. Another interesting observa- 
tion is the shell structure exhibited by all finite-size (fs) 
results. In order to visualize both the fixed-node error 
and the finite-size effects, we have plotted in Fig. Q] the 
following calculations: the static density-response func- 
tion of (i) an infinitely large unpolarized non-interac ting 
free electron gas, the well-known Lindhard function [ll| 
Xo (dashed line), (ii) Xo S of a non-interacting unpolar- 
ized system of 114 and 4218 electrons (solid lines), (iii) 
an infinitely large unpolarized interacting free electron 
gas in the random-phase approximation (RPA), xrpa 
(dotted line), (iv) a finite unpolarized system of 114 non- 
interacting electrons within the fixed-node approxima- 
tion giving Xo (top dot-dashed line), (v) a finite un- 
polarized system of 114 interacting electrons within our 
fixed- node DMC scheme (middle dot-dashed line), and 
(vi) a finite unpolarized system of 114 interacting elec- 
trons in the fn RPA x^rpa (dot-dashed line at the bottom, 
see below for details). We see that the finite-size shell 
structure is not negligible even for a system of 4218 elec- 
trons. Nevertheless, the exact non-interacting density- 
response function nicely reproduces the well-known Lind- 
hard function especially for k/kf > 2. 

C. Discussion 

Since the fn Lindhard function Xo™ (top dashed line in 
Fig- HI is smaller than the real Lindhard function xo (dot- 
ted line at the top of Fig. [2]) the fn interacting x^™ is also 
too small (lower solid line in Fig. [2] compared to the dots 
showing the relevant data of Ref. [2]). However, assuming 
that the effect of the fixed-node nature of the calculation 
is the same for the interacting and non-interacting case, 
we should still be able to extract meaningful correlation 
data and reverse the effects of the fn approximation [8j . 
Let us start with the definition for the exchange correla- 
tion (xc) kernel f xc and the local field factor G 

- f xc (k) = v c {k)G(k) = J- l — 

X( k ) XRPA(k) 
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FIG. 2: The density response at r 3 = 5 for an unpolarized 
114 particle homogeneous electron gas. The dotted lines show 
the Lindhard function xo (top) and \Rpa (bottom). The 
dashed lines are Xo (top) and Xrpa (bottom). Finally, the 
solid lines show the fixed-node x^" (bottom) and extrapolated 
interacting x- The dots correspond to the results in Ref. 0|. 
Our data reaches to smaller k values as our system is larger 
(114 vs. up to 66 electrons). 



FIG. 3: The ratio 



Xii 



at r s 



5 for TV = 66 and N 



X(k) Xo(k) 



v c (k) 



(10) 



where Vc(k) = |f . In practise, we do not have access to 
x(k) but only to its finite-size equivalent in the fixed-node 
approximation. Moroni and coworkers who include 
the nodal variation at a Kohn-Sham level argue that that 
while the density response contains finite size effects, f xc 
is less afflicted by these. Hence, they extract f xc and add 
that back on to the non fixed-node infinite-cell Lindhard 
function to correct for finite size effects, thus eliminat- 
ing the shell structure. There is no reason to expect 
the nodal variation of the KS nodes to correctly describe 
the nodal variation of the fully interacting system with 
respect to \: The KS nodes and the true many-body 
nodes are unrelated. Furthermore, different QMC sys- 
tems at different numbers of electrons N also correspond 
to distinct nodes, but the data for f xc (e.g. Ref. [2] or 
Fig. [7J for different values of N is mutually compatible. 
Finally, the effect of the nodal variation on xo seems uni- 
versal, i.e. independent of N except for shell effects (see 
Fig. [3]). It therefore seems reasonable to assume that f xc 
is independent of nodal effects and can thus be used to 
correct for wrong or absent nodal variation. Hence, by 
using the fixed-node quantities in Eq. (|10p . implicitly 
defining a fixed- node Xrpa> we can derive a fn and 
Qfn (Figures U]and[H]). These are remarkably similar to 
the data in Ref. B. I n fact, //" even has a slight dip as 
suggested in Ref. [12| which however is not really visible 
in Ref. [2]. This is encouraging and indeed we can use 
our data for G? n in conjunction with the real Xrpa in 
Eq. (fTU)) to estimate the non-fn interacting x- The result 



114 electrons. Due to scaling this graph is independent of 
the value of r s , except for noise. Note the pronounced shell 
structure for k < 2k f. 



can be seen in Fig. [2] All the fn quantities are too small 
compared to their non-fn counterparts. However our ex- 
trapolated date (solid line at the top) is very close the 
extrapolated data of Ref. Q (dots). 

For completeness sake Fig. [6] shows details of the 
extrapolated x at r s — 2,5, and 10. Also, in Fig. [7J we 
show a direct comparison between our results and Ref. Q 
where it is possible, i.e. at N — 66 electrons in addition 
to our results for N = 114 electrons confirming that all 
our data is compatible with Ref. @- Except for noise 
there is no significant difference between data at different 
N, corroborating the assumption that f xc is independent 
of finite-size effects. 

In general, our results nicely follow the data in Ref. , 
who take into account the change in the nodes at a Kohn- 
Sham level, whereas our calculations do not take into ac- 
count any nodal effect on xc quantities. The fact that 
the methods yield consistent results for f xc suggest that 
assuming f xc to be free from nodal effects is justified and 
that in either case the resulting data is an accurate de- 
scription of systems with the full interacting nodal vari- 
ation. 



IV. 



CONCLUDING REMARKS 



We have generalized Ref. [6j to the second derivative 
of the energy. This yields a novel method and an efficient 
algorithm to calculate the static response function within 
DMC. Our algorithm permits the computation of a large 
number of diagonal and off-diagonal terms in a single 
DMC run without the need for numerical derivatives or 
re-optimization. Noise can be efficiently controlled by in- 
creasing the number of DMC walkers and we have found 
that we can use large DMC time steps without introduc- 
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FIG. 4: /i c n at r s = 2, 5, 10. The value at k = 0.4189fc F 
clearly is an outlier. Also, the noise increases as k grows. 
Interestingly, there is a slight dip in //" for k/kF < 2 as 
demonstrated in Ref. [Q. The dots correspond to the data 
in Ref. 0. 




FIG. 5: As figure [4j but showing the local field factor. 



ing a bias, potentially speeding up calculations greatly. 
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